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For precision studies with QCD observables at colliders, higher order perturbative cor- 
rections are often mandatory. For exclusive observables, like jet cross sections or differ- 
ential distributions, these corrections were until recently only known to next-to-leading 
order (NLO) in perturbation theory. Owing to many new technical developments, first 
next-to-next-to-leading order (NNLO) QCD calculations are now becoming available. 
We review the recent progress in this field. 

1 Introduction 

At present-day colliders, a number of benchmark reactions is measured to a very high experi- 
mental accuracy. These reactions allow for a very precise determination of the parameters of 
the Standard Model, and may reveal minute deviations, hinting towards new physics effects. 
For many of these benchmark reactions, one faces the problem that the observables [1] are 
inherently exclusive (like jet cross sections), or differential in several variables (if kinemat- 
ical cuts are involved). While fully inclusive quantities arc often known to third order in 
perturbation theory, the perturbative expression for those exclusive observables is much less 
well known. Consequently, the precision of many benchmark reactions is limited not by the 
quality of the experimental data, but by the error on the theoretical (next-to-leading order, 
NLO) calculations used for the extraction of the Standard Model parameters. To improve 
upon this situation, an extension of the theoretical calculations to next-to-next-to-leading 
order (NNLO) is therefore mandatory. 

2 Ingredients to NNLO calculations 

At next-to-next-to-lcading order, three types of processes contribute to an observable with 
m hard particles in the final state: 

1. the two-loop corrections to the m-particle process, where all particles are hard. 

2. the one- loop corrections to the (m+ l)-particle process, where one particle can become 
unresolved (collincar or soft). 

3. the tree-level (m+2)-particle process, where up to two particles can become unresolved. 

Each contribution contains infrared singularities, which cancel only in the sum. Since the 
definition of the observable final state (jet algorithm or kincmatical cuts) acts differently 
on each of these processes, it is mandatory to compute all of them individually. This is 
in variance with calculations of fully inclusive quantities (total cross sections, sum rules), 
where all contributions can be added before evaluation using the optical theorem. 

QCD infrared factorisation predicts that the behaviour of cross sections in soft or collinear 
limits is process- independent. At NNLO, the corresponding universal tree- level double un- 
resolved [2] and one-loop single unresolved [3] factors are known. Based on this universal 
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behaviour, one aims to develop process-independent techniques for NNLO calculations of 
exclusive observables. 

The first calculations combining elements of this type are the derivations of the NNLO 
corrections to the deep inelastic structure functions [4], and to the coefficient functions for 
vector boson [5] and Higgs boson production [6,7] at hadron colliders. Although these ob- 
servables are fully inclusive, one still has to consider the individual parton-level contributions 
separately, since they include collincar radiation associated with the incoming partons. 

NNLO calculations for hadron colliders also require parton distributions accurate to this 
order. Following the derivation of the three- loop splitting functions [8], global NNLO fits [9] 
are now becoming available. These are however still somewhat limited by the fact that not 
all observables included in the fit are known to NNLO. 

3 Techniques and applications 

The derivation of the individual ingredients to NNLO calculations, and their combination 
into a parton-level event generator, are posing a variety of computational challenges. Many 
new techniques have been developed in this context. 

3.1 Sector decomposition 

In computing the various contributions to perturbative corrections at NNLO and beyond, 
one frequently encounters the problem of overlapping singularities. These appear if one 
particular term develops a singular behaviour in more than one region of phase space. The 
technique of sector decomposition offers an elegant solution to this problem. Originally 
proposed as a tool in formal proofs of renormalisability [10], and only used occasionally 
afterwards [11], this technique was fully formulated first for multi-loop integrals [12] and 
subsequently extended to phase space integrations [13, 14]. 

Starting from a parameter representation of the (virtual or real) integral under consid- 
eration, the space of integrations is decomposed iteratively into non-overlapping sectors. 
At the end of this iteration each sector contains only a single type of singularity. In each 
sector, after expanding all regulators in distributions, the Laurent series of the integral is 
well-defined, and its coefficients can be obtained by numerical integration. 

This technique has been used in the calculation of virtual two-loop and three-loop multi- 
leg integrals [12,15], often ahead of or in timely coincidence with analytical results. For loop 
integrals involving mass thresholds inside the physical region, sector decomposition must be 
combined with contour deformation to avoid pinch singularities [16,17]. 

Concerning the calculation of real radiation corrections at NNLO, the first-ever results for 
fully differential observables were obtained using sector decomposition for e + e~ — > 2 jets [18] , 
pp^ H + X [19], muon decay [20] and pp -> V + X [21]. 

3.2 Integral reduction 

Within dimensional regularisation, the large number of different integrals appearing in multi- 
loop calculations can be reduced to a small number of so-called master integrals by using 
integration-by-parts (IBP) identities [22]. These identities exploit the fact that the integral 
over the total derivative of any of the loop momenta vanishes in dimensional regularisation. 
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For integrals involving more than two external legs, another class of identities exists 
due to Lorentz invariance. These Lorentz invariance identities [23] rely on the fact that 
an infinitesimal Lorentz transformation commutes with the loop integrations, thus relating 
different integrals. The common origin of IBP and LI identities is the Poincare invariance 
of loop integrals within dimensional regularisation. 

In principle, the identities for all loop integrals of a given topology can be solved in 
a closed symbolic form [22]. In practise, this solution is often overly complex, and can 
not be found in an automated manner. Instead, one pursues another approach, the so- 
called Laporta algorithm [24], which solves the system of identities for a given topology by 
assigning a weight to each integral according to its complexity. In turn, the very large system 
of interconnected identities is solved by elimination, thereby reducing each integral in the 
system to master integrals. Several computer implementations of the Laporta algorithm are 
available [23-25]. 

Originally formulated only for loop integrals, these reduction techniques can also be ap- 
plied to multi-particle phase space integrals by expressing on-shcll conditions and kinematic 
constraints in the form of generalised propagators [7, 14] . 

3.3 Mellin-Barnes integration 

The techniques described in the previous section allow to reduce the hundreds to thousands 
of different integrals appearing in an actual calculation to a small set of master integrals. 
The reduction equations do not yield any information on the master integrals, which have 
to be determined using some other technique. Unlike in the one-loop case, where a direct 
integration using Feynman parameters is usually sufficient, master integrals at two or more 
loops pose a considerable computational challenge. 

A commonly used analytical approach is to replace the product of loop propagators by 
multiple Mellin-Barnes integrations [26]. This replacement can be performed in an auto- 
mated manner using computer algebra [27]. After integration over the loop momenta in 
dimensional regularisation, the Mellin-Barnes integrals are normally not yet well-defined 
around d = 4, where the Laurent expansion is desired. To arrive at this limit, one must first 
perform analytic continuations, which can again be done using computer algebra [28,29], 
thereby transforming each single Mellin-Barnes integral into a sum of integrals and residues. 
After analytic continuation, the Laurent coefficients can be determined by carrying out the 
Mellin-Barnes integrals numerically [28]. Analytical results for Mellin-Barnes integrals can 
often be obtained (especially if they can be expressed as multiple harmonic sums [30-32]), 
and systematic methods for them are under development. 

Using Mellin-Barnes integration techniques, analytical results were obtained for two-loop 
four-point functions with massless [33-35] and massive [36] internal propagators. 

3.4 Differential equations 

A method for the analytic computation of master integrals avoiding the explicit integration 
over the loop momenta is to derive differential equations in internal propagator masses or 
in external momenta for the master integral, and to solve these with appropriate boundary 
conditions. This method has first been suggested [37] to relate loop integrals with internal 
masses to massless loop integrals. 

It has been worked out detail and generalised to differential equations in external mo- 
menta in [23,38]. Differentiation of a master integral with respect to an external invariant 
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or to a mass yields a combination of integrals of the same topology. Using the integral 
reduction techniques described above, these can be reexpressed by the master integral it- 
self, plus simpler integrals. As a result, one obtains an inhomogeneous differential equation 
for the master integral. Solving this equation and matching the solution onto an appro- 
priate boundary condition (obtained in a special kinematical point) then yields the desired 
master integral, very often in a closed form containing hypergeometric functions and their 
generalisations. The Laurent expansion of the integral then amounts to expansion of these 
functions [39] . A detailed review of the method can be found in [40] . 

Using the differential equation technique, master integrals were derived for massless two- 
loop four-point functions [41-43], for massive two-loop vertex functions [44,45] and for parts 
of the master integrals required for Bhabha scattering [46] . 

3.5 Virtual corrections and infrared structure 

Using the IBP and LI equations for the reduction to master integrals, which were then 
computed with various of the above-mentioned techniques, results for a variety of two-loop 
corrections were obtained for 1 — > 2, 1 — ► 3 and 2 — ► 2 reactions. These include all massless 
parton-parton scattering amplitudes [47], processes yielding two-photon final states [48], 
light-by-light scattering [49] and vector boson decay into three massless partons [50,51] and 
its crossings [52]. For amplitudes involving external masses, all two- loop form factors [53] of 
a heavy quark were derived, and large parts of the virtual two-loop corrections to Bhabha 
scattering [54] are completed. 

These results are usually expressed in terms of harmonic polylogarithms [55] , which are 
a generalisation of the well-established Nielsen's polylogarithms [56]. If several scales are 
involved in a process, the set of functions needs to be further extended to include multi- 
dimensional harmonic polylogarithms [50,57,58]. 

After ultraviolet renormalisation, these two-loop amplitudes still contain poles of infrared 
origin. This infrared pole structure is universal for massless amplitudes, and can be predicted 
from resummation formulae [59]. Exploiting the fact that a particle mass can also act as 
an infrared regulator, the universality of infrared singularities can be extended to massive 
amplitudes [60] , where not only the divergent terms but also logarithmically enhanced terms 
are universal. Knowing the corresponding massless scattering amplitudes, it is therefore 
possible to construct massive amplitudes up to corrections of order m 2 /s, which was recently 
accomplished for heavy quark production [61] and Bhabha scattering [62,63]. 

3.6 Subtraction methods 

To build exclusive final states at a given order, a jet algorithm or event shape definition has 
to be applied separately to each partonic channel contributing at this order and all partonic 
channels have to be summed. However, each partonic channel contains infrared singularities 
which, after summation, cancel among each other. Consequently, these infrared singularities 
have to be extracted before the jet algorithm can be applied. While explicit infrared sin- 
gularities from purely virtual contributions are obtained immediately after integration over 
the loop momenta, their extraction is more involved for real radiation. The singularities 
associated with the real emission of soft and/or collinear partons in the final state become 
only explicit after integrating the real radiation matrix elements over the appropriate phase 
space. 
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In the sector decomposition method described above, the full matrix elements are inte- 
grated numerically, after expansion as a Laurent series (which amounts to subtraction of 
residues in each sector). A different approach is pursued by subtraction methods, which 
extract infrared singularities of the real radiation contributions using infrared subtraction 
terms. These terms are constructed such that they approximate the full real radiation matrix 
elements in all singular limits while still being integrable analytically. 

Several methods for constructing NLO subtraction terms systematically were proposed 
in the literature [64-68] . For some of these methods, extension to NNLO was discussed [69] 
and worked out for special cases [70]. Up to now, the only method worked out in full detail 
to NNLO is antenna subtraction [71]. 

The basic idea of the antenna subtraction approach is to construct the subtraction terms 
from antenna functions. Each antenna function encapsulates all singular limits due to the 
emission of unresolved partons between two colour-connected hard partons. All antenna 
functions can be derived systematically from matrix elements [72] for physical processes. 
The antenna subtraction method was used recently in the derivation of the NNLO QCD 
corrections to e + e~ — ► 3 jets and related event shapes [73]. 



4 Results and outlook 

The computational challenges of NNLO QCD calculations required for precision phenomenol- 
ogy have led to the development of a variety of new technical methods. Using these methods, 
many ingredients to NNLO QCD calculations were assembled in recent times. Many core 
results were derived more than once using independent, different methods. 

These ingredients are now assembled into complete NNLO QCD calculations, which are 
usually carried out in the form of a parton-level event generator. In such a program, the 
full kincmatical information is available for each event, and cross sections are obtained by 
adding all events relevant to the observable under consideration. This setup allows a great 
flexibility in the implementation of experimental cuts, and in the simultaneous evaluation 
of numerous kincmatical distributions. 

Up to now, fully exclusive NNLO QCD calculations were performed for e + e~ — > 2j [18, 
72,74], e + e~ — > 3j [73,75], the forward-backward asymmetry in e + e~ annihilation [76], 
Higgs production [19, 70] and vector boson production [21] at hadron colliders, as well as 
muon decay in QED [20] . 

Virtual two-loop corrections are available for a number of further 2 — > 2 reactions, includ- 
ing jet production and vector-boson-plus-jet production at hadron colliders, jet production 
in deep inelastic scattering, and heavy quark production observables. The full calculation 
of these observables may require further improvements to existing tools, but further results 
of NNLO calculations are clearly within reach in the near future. 
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